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Abstract —Localization of radio frequency sources over mul¬ 
tipath channels is a difficult problem arising in applications 
such as outdoor or indoor gelocation. Common approaches that 
combine ad-hoc methods for multipath mitigation with indirect 
localization relying on intermediary parameters such as time-of- 
arrivals, time difference of arrivals or received signal strengths, 
provide limited performance. This work models the localization of 
known waveforms over unknown multipath channels in a sparse 
framework, and develops a direct approach in which multiple 
sources are localized jointly, directly from observations obtained 
at distributed sources. The proposed approach exploits channel 
properties that enable to distinguish line-of-sight (LOS) from 
non-LOS signal paths. Theoretical guarantees are established 
for correct recovery of the sources’ locations by atomic norm 
minimization. A second-order cone-based algorithm is developed 
to produce the optimal atomic decomposition, and it is shown to 
produce high accuracy location estimates over complex scenes, 
in which sources are subject to diverse multipath conditions, 
including lack of LOS. 


I. Introduction 

Traditional time-of-arrival (TOA)-based localization is ac¬ 
complished through a two-step process. In the first step, 
sensors estimate TOA’s from all incoming signals; in the 
second step, such estimates are transmitted to a central node, 
that subsequently estimates the location of each source by 
multilateration [1], We refer to these localization techniques 
as indirect localization. In a multipath environment, each 
sensor receives, in addition to a line-of-sight (LOS) signal, 
multiple (possibly overlapping) replicas due to non-line-of- 
sight (NLOS) paths. Due to these multiple arrivals, it is, in 
general, more challenging to obtain accurate TOA estimates 
of the LOS components at the sensors. Matched filtering is a 
method for time delay estimation. However, its performance 
degrades greatly in the presence of multipath whose delay is 
of the same order than the inverse of the bandwidth of the 
signal [2]. Moreover, in the case of blockage of the LOS path, 
the TOA of the first arrival does not correspond to a LOS 
component anymore, and will corrupt localization. In such a 
case, it is customary to apply techniques, like the one in [3], 
to mitigate NLOS channel biasing of the geolocation estimate. 

A better approach than indirect localization is to infer 
the source locations directly from the signal measurements 
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without estimating any parameters such as propagation delays. 
The concept of direct localization was first introduced by 
Wax and Kailath [4, 5] in the 70’s. However, it is in the 
last decade that Weiss et al. have further investigated and 
proposed actually efficient techniques [6-9]. In the absence of 
multipath, the state-of-the-art is Direct Position Determination 

(DPD) [6] which outperforms standard indirect localization, 
particularly at low signal-to-noise ratio (SNR), because it takes 
into account the fact that signals arriving at different sensors 
are emitted from the same location. The literature on direct 
localization in the presence of multipath is scarce. In [10] a 
maximum likelihood (ML) estimator has been developed for 
the location of a single source assuming a fixed and known 
number of multipath, but without providing an efficient way to 
compute the estimator. In [9], a Direct Positioning Estimation 

(DPE) technique is proposed for operating in dense multipath 
environments, but requires knowledge of the power delay 
profile and is limited to localization of a single source. 

A requirement of direct localization is that the signals, 
or a function of them, are sent to a fusion center which 
estimates the source’s locations. Thus direct techniques are 
best suited for centralized networks. An example of this are 
Cloud Radio Access Networks (C-RAN) [11, 12], C-RAN is 
a novel architecture for wireless cellular systems whereby the 
base stations relay the received signals to a central unit which 
performs all the baseband processing. Cellular systems are 
required to be location-aware, that is they must be able to 
estimate the locations of the user equipments (UE) for appli¬ 
cations such as security, disaster response, emergency relief 
and surveillance in GPS-denied environment [13]. In addition, 
in the USA, it is required by the Federal Communications 
Commission (FCC) that by 2021 the wireless service providers 
must locate UE’s initiating an Emergency 911 (E-911) call 
with an accuracy of 50 meters for 80% of all wireless E-911 
calls [14], In uplink localization, the base stations perform 
time measurements of the received signals emitted by the UE’s 
in order to infer their positions. Thus, a high accuracy direct 
localization technique designed for multipath channels, such 
as the one proposed in this work, may enhance the localization 
accuracy of current existing cellular networks by utilizing the 
C-RAN infrastructure. Moreover, it exists other applications 
that may benefit from high accuracy TOA-based geolocation 
such as in WLAN and WPAN networks. For instance, the 
setup of [15] uses radios with the IEEE 802.15 (WPAN) 
standard to localize devices. The setups proposed in [16, 17] 
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employ TOA-based localization for localizing 802.11 devices 
(WLAN). In [18] it is proposed a hybrid RSS(received signal 
strength)-TOA based localization algorithm that works with 
802.11 and 802.15 technologies. Other TOA-based localization 
applications are in the radio frequency identification (RFID) 
field [19]. 

In this paper, we present a TOA-based direct localization 
technique for multiple sources in multipath environments 
(DLM) assuming known waveforms and no prior knowledge of 
the statistics of the NLOS paths. Preliminary results of DLM 
were presented in [20]. Without some prior knowledge on the 
multipath, NLOS components carry no information, and the 
best performance is obtained by using only LOS components 
[13]. We propose an innovative approach, based on ideas of 
compressive sensing and atomic norm minimization [ 21 ], for 
jointly estimating the sources’ locations using as inputs the 
signals received at all sensors. Numerical evidence shows 
that DLM has higher accuracy than indirect techniques, and 
that it works well in a wide range of multipath scenarios, 
including sensors with blocked LOS. Moreover DLM requires 
no channel state information. 

In Section II, we introduce the signal model. Section III 
briefly introduces our proposed technique. Sections IV and 
V provide in-depth explanations on the different parts of our 
technique. Section VII compares DLM to previous existing 
techniques. Finally, Section VIII reports our conclusions. 

II. Signal Model 

Consider a network composed of L sensors and Q sources 
located in a plane. The location of the q- th source is defined by 
two coordinates stacked in a vector p, r All sources share the 
same bandwidth B, and transmit their own signals {s q (t)Y q= i- 
The number of sources Q and their waveforms are known. 
The observation time is T, assumed to be shorter than the 
time coherence of the channel, therefore, the channel is time- 
invariant. The complex-valued baseband signal at the /-th 
sensor is 

n(t) = rl os (t) + rf LOS (t) + w l (t) 0 <t<T, (1) 

where wi(t) is circularly-symmetric complex white Gaussian 
noise with known variance E \wi(t )\ 2 = a 2 ,. The term r\ os (t) 
is the sum of all LOS components: 

Q 

r\- os (t) = ^2a q is q {t-Ti(p q )), (2) 

9=1 

where a q i is an unknown complex scalar representing the 
signal strength and phase of the LOS path between the q- 
th source and /-th sensor, and 17 (p) is the delay of a signal 
originating at p and reaching the /-th sensor: 

t i(p) = Up - PtII 2 / c - (3) 

In (3), pj is the location of the /-th sensor, c is the speed of 
light and || • H 2 denotes the standard Euclidean norm. The term 
rf LOS (t) in (1) aggregates all NLOS arrivals: 

Q M q , 

= EE (* - r<r>) . (4) 

q— 1 m—1 


where M q i denotes the unknown number of NLOS paths 
between the q- th source and the /-th sensor, a^ is an 
unknown complex scalar representing the amplitude of the 
rn-th NLOS path between the r/-th source and /-th sensor, 
and rjj " 1 is the delay of the NLOS component. The received 
signal (1) is sampled at a frequency f s satisfying the Nyquist 
sampling criterion: f s > 2 B, where B is the bandwidth of 
r(t). Each sensors collects N time samples at each observation 
time. By stacking the N acquired samples, the received signal 
r 1 = [r;(0),... ,ri((N — 1 )/f s )] T at the /-th sensor can be 
written in the following vector form 

Q Q M q i 

n = £ «9* S 9 MPq)) + £ H a S i)s 9 ( T g ( / m) ) + Wi > (5) 

q— 1 q =1 m—1 

where s q (r) is the vector of the N received samples from the 
g-th source waveform with delay r: 

S 9 ( T ) = [ S 9(°~ T ) ••• S 9 i( N - !)//« - T )] 1 ' ( 6 ) 

Since all sensors acquire the same number of samples, the 
samples may be stacked in an N x L matrix 

R = [ri • • • r L ] = 

Q 

= Kl S 9 (ti(P<j)) ••• U qL S q (T L {p q ))\ + 

9=1 

Q L M ql 

+ EEE“f )s 9 ( r 9 ( r ) ) V f + W > (7) 

q— 1 l—l m—1 

where the rows and columns index time instants and sensors, 
respectively, and V; is an all-zeros vector except for the Z- 
th entry which is one. The LOS and NLOS components 
are parametrized by the first and second summands in (7), 
respectively. In the rest of the paper, we will switch between 
the notations in (5) and (7) depending on whether we are 
interested in the signal of one sensor only or of all sensors. 

III. Proposed Localization Technique 

In order to develop a localization technique, it is first 
necessary to understand what parameters of the received 
signals depend on the sources locations. In the signal model 
introduced in the previous section, the propagation delays of 
the NLOS components (4) were assumed to be unknown and 
arbitrary, because of the lack of prior statistical knowledge 
of the channel. Thus, information on the sources locations 
is carried only by the LOS components (2). This claim is 
supported by the analysis in [13], which showed that the CRB 
increases when NLOS components are present. Consequently, 
without a priori knowledge, the optimal strategy is to reject 
NLOS components as much as possible, and rely on the LOS 
components to infer the sources’ locations. 

With indirect techniques, first, the TOA’s of the LOS com¬ 
ponents are estimated, and then used to localize the sources by 
multilateration. However, indirect techniques are suboptimal 
because they estimate the TOA of the first path at each sensor 
independently, instead of taking into account that all LOS 
components originate from a single source location. In this 



3 


section, we propose a direct localization technique that relies 
on the fact that all LOS components associated with a source 
must originate from the same location. Under the Gaussian 
assumption, the maximum likelihood estimator (MLE) is the 
solution to the following fitting problem 


mm ) 
Pl.---.PQ /. - J 

an,...,a L Q l—l 


.(l) 


LQ 

( m LQ ) 
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J 2 a ql S q ( T l(P q )) ~ 

9=1 

Q M q , 
q— 1 m— 1 
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( 8 ) 


subject to > r/(p g ), for all q, l and m. The parameters of 
interest are the source locations {Pq}^ =1 , while the rest act as 
nuisance parameters. Besides the fact that it is an enormous 
challenge to find an efficient technique for minimizing this 
objective function, the ML criterion does not even lead to a 
satisfactory solution. The reason is that M q i, for all Z and 
q , are hyperparameters that control the number of NLOS 
paths in our model. It is known that increasing the values 
of hyperparameters always leads to a better fitting error [22], 
and in our case, it would lead to the erroneous conclusion 
that there are a very large number of NLOS arrivals. Instead, 
we assume that the number of NLOS arrivals and the number 
of sources is low with respect to the number of observations. 
This assumption enables the formulation of a feasible solution 
to the ML multipath estimation problem by means of a sparse 
recovery technique. 

In order to obtain a high-precision localization technique, 
there are two properties of the signal paths that need to be 
exploited. These properties allow to distinguish LOS from 
NLOS components. The first one is that NLOS components 
arrive with a longer delay than LOS components, and the 
second property is that all LOS paths originate from the same 
location. Our technique is divided into two stages, which are 
explained in the following two sections. In the first stage, 
NLOS components are canceled out from the received signals 
by exploiting the fact that LOS components must arrive first. 
This processing can be done locally at each sensor. In the 
second stage, the cleaned version of the received signals are 
sent to a fusion center that finds the sources’ location. It is in 
this stage that the source locations are estimated by exploiting 
the fact that LOS components must originate from the same 
location, whereas NLOS components may be local to the 
sensors. 


IV. Stage 1: Deconvolution 

In this stage, the multipath channel is deconvolved, or 
equivalently, the propagation delays of different paths are 
estimated, and the multipath contributions are removed from 
the received signals. Our technique of choice for deconvolution 
is the sparsity-based delay estimation technique proposed by 
Fuchs [23] because of its high accuracy and because it uses 
only a single snapshot of data as in our case. Other high 
accuracy time delay estimation methods, like MUSIC [24], are 
not applicable here because they require multiple uncorrelated 
data snapshots. Let r max be the largest possible propagation 


delay, then in the Fuchs’ technique, the continuous set of all 
possible propagation delays [0, r max ] is discretized forming a 
grid of delays 

V = {0,T res , • • • ,r max } , (9) 


where parameter T res denotes the resolution of the grid. Define 
the dictionary matrix stacking the received signal waveforms 
for all possible (discrete) delays (9): 

A = [s(0) ••• s(r max )] . (10) 


Then, the propagation delays of all paths from a single source 
reaching the Z-th sensor are estimated by solving the following 
Lasso problem of the form 

mmA Hxlh + ||rj - Ax^ , (11) 

where A is a regularization parameter, || • ||i is the £\- 
norm of a vector, and r/ is the received signal defined in 
(5). Solving this convex optimization problems, results in a 
sparse vector x whose non-zero entries indicate the estimated 
delays. More precisely, if the cZ-th entry of x is different 
than zero, then a path has been detected with propagation 
delay (d — l)T res . After estimating the propagation delays, 
Fuchs uses a maximum description length (MDL) criterion to 
filter out false detections. For more details on this technique 
see [23], and for a better understanding on the mathematics 
behind the Lasso problem see [25], In [23], the time delay 
estimation technique was designed for real-valued signals, and 
assuming only a single emitting source. Here, we generalize 
such approach to complex valued signals by simply allowing 
the variables and parameters in (11) to be complex. We also 
generalize it to multiple sources by expanding the columns of 
the dictionary (10) to the waveforms of all sources: 


A = [si (0) • • • si (r max ) • • • sq (0) • • • s Q (r max )] . 

( 12 ) 

It is possible to use other delay estimation techniques. Ob¬ 
viously, the more accurate the delay estimation technique, 
the better performance would be expected from this NLOS 
interference mitigation. Contrary to indirect localization tech¬ 
niques, the goal here is not to precisely estimate the prop¬ 
agation delays of the first paths, but rather to estimate the 
propagation delays of all subsequent arrivals, and cancel them 
out. 

Let f ql ,..., t^" 1 be the estimated propagation delays from 
source q to sensor Z, then their amplitudes may be estimated 
by solving a linear least squares fit 


{<■} 


= argmin 

K«} 


Q Pql 
9=1 P= 1 


a p s 


( f j) 


(13) 


Assuming the estimated propagation delays are ordered in 


ascending order f ql 


< 


< T, 


P,l 


ql 


then all arrivals, except 


the first, can be canceled out from the received signals 


Q Pql 


~ r l= V l~Y,Y, 5 9* S 9 i f ql) 


(14) 


9= 1 P =2 


Ideally, all NLOS arrivals would be perfectly detected and 
their propagation delays estimated, in which case we could 
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continue with a direct localization technique designed for 
absent multipath. However, (14) is not guaranteed to cancel 
all NLOS components for two reasons. First, if the LOS 
path between a source and sensor is blocked, then the first 
arrival corresponds to a NLOS paths, in which case it is not 
removed. Also, it is possible that the chosen delay estimation 
technique misses some arrivals or detects some false ones, thus 
failing to remove some NLOS components or adding some 
extra components, respectively. In short, this stage is essential 
as it reduces the multipath, but does not necessarily remove 
it completely. In the next section, we present a localization 
technique designed to work in the presence of the residual 
multipath as well as blocked paths. 

V. Stage 2: Localization 

This stage seeks to estimate the sources locations using 
the signals {r i}j" =1 output by Stage 1. As explained in the 
preceding section, such signals include LOS and also NLOS 
components, therefore, the signal model introduced in (5) for 
r i is also valid for r/. Obviously, since f; and r; are different, 
so are the values of the parameters appearing in (5) that 
characterize them. From here on, to keep the notation in check, 
we abuse the notation by writing r; instead of f/. However, 
always bear in mind that the observations in this stage are the 
signals output by Stage 1. 

To compute the MLE (8), it is required that the number of 
LOS and NLOS paths be known, otherwise the minimization 
(8) tends towards a nonsensical solution with a very large 
number of paths, a problem known as noise overfitting [22]. In 
this section, it is first assumed that the number of sensors that 
receive a LOS path from the g-th source, say S q , is known. It 
will be shown later in Section V-C that such information is not 
really needed. Nevertheless, even if {,S' g }^ =1 are known, but 
since the number of NLOS paths is not, a pure MLE approach 
is still not feasible. To bypass this issue, we will rely on the 
fact that the number of sources and NLOS paths is relatively 
small with respect to the number of observations. 

Define a LOS atom as the N x L matrix of measurements 
of LOS paths of a signal s q (t) emitted from location p and 
received at the L sensors, i.e., 

L<? (b,p) = [b(l)s q (n(p)) ••• b(L)s q (tl(p))] (15) 

where b = [6(1) • • • b{L)] T are the complex amplitudes of the 
LOS components. It is important to normalize b as it will be 
discussed shortly. Hence, ||b ||2 is constrained to a given value 
that we will denote u q , i.e., ||b ||2 = u q . Define a NLOS atom 
as the N x L matrix of measurements due to a single NLOS 
path from the q- th source to the /-th sensor 

( r ) = e ^ s <?( T ) v r ( 16 ) 

where the phase and delay are <j> and r, respectively, and v/ 
is a unit vector, with the unit entry indexed by l. Note that 
the dependence of ~N q i (r) with respect to f is omitted in the 
notation for simplicity reasons. Let R denote the matrix of 
received signals (7) in the absence of noise. Then, R may be 
expressed as a positive linear combination of given atoms 

c ( fc) A ( fc ), A (k) £A (17) 

k 


where c <k ' 1 > 0 for all /c, and A is the set of all atoms (or 
atomic set). The atomic set includes all different LOS and 
NLOS atoms, 

A = -4los U ,4nlos (18) 


where Alos 


-4los= U {L g (b )P ):beC L ,pe5cM 2 ,||b|| 2 =u g } 
q =1 

(19) 

and LIxlos 


Q L 

-4nlos = U U { N «' ( t ) : 0 < <£ < 2 tt, r G [0,T max ]|. 

q=l 1=1 

( 20 ) 

Here, S denotes the search area of the sources and r max the 
maximum possible delay in the system. Notice that the set of 
LOS atoms and the set of NLOS atoms are infinite in the sense 
that p and r are continuous variables within their domains. 
Thus this framework is inherently different in that the discrete 
dictionaries used in traditional compressive sensing. 

Since the atomic sets are infinite, determining the coeffi¬ 
cients c <k> from measurements R is a highly undetermined 
problem. This problem is solved by seeking a sparse or simple 
solution in some sense to the coefficients c <k> . As motivated 
in [21], this can be accomplished with the help of the concept 
of the atomic norm. More precisely, the atomic norm || • ||_4 
induced by the set A is defined as 



3 (fc) :R=^c ( ' !) A(‘),A (fe) gA\. 


( 21 ) 


An atomic decomposition of R is any set of coefficients {c^ fc l} 
for given atoms {AW} such that R = f2k c (fc)A( fc ). The 
cost of an atomic decomposition is defined as the sum of 
its positive coefficients: atomic decomposition is 

optimal if its cost achieves |jR||_ 4 , or equivalently, if its cost 
is the smallest among all atomic decompositions. Sparsity is 
imposed here in the sense that we assume that the coefficients 
c!' k> for which the atomic decomposition is optimal (i.e., lowest 
cost) are associated with the true solution of locations and 
time delays. This sparsity condition resolves the undetermined 
nature of (17). In practice, in the presence of noise, we seek 
the optimal atomic decomposition that approximately matches 
the received signals. Precisely, in [21] it is suggested that the 
noiseless signals R may be estimated by minimizing 


min 

R 1 

R 

L 

2 

(22a) 

S.t. 

R 

- R 

< e, 

F 

(22b) 


where || • \\p is the Frobenius norm, i.e., ||M||/ = 

\J^2i j \M(i,j)\ 2 , f° r an y matrix M, where M(i,j) is the 
entry at the Lth row and j-th column. Roughly speaking, 
minimizing the atomic norm (22a) enforces sparsity, while 
constraint (22b) sets a bound on the mismatch between the 
noisy signals and the estimated signals. In fact, the left hand 
side of (22b) is an ML-like cost function (8), hence, parameter 
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e may be regarded as an educated guess of the ML cost. The 
optimum solution to problem (22), say R*, may be regarded 
as an estimate of the received signals in the absence of noise. 
However, notice that solving such problem produces R* only 
and not its optimal atomic decomposition. Thus, in general, in 
order to recover the optimal atomic decomposition, first, the 
optimum R* to problem (22) is computed, and second, the 
optimal atomic decomposition of R* is found. 

The atomic decomposition of R* may be expressed 

R* = EX>fL, ( b W)+E EX>>* {<?) 

q—1 k—1 q= 1 /—1 k—1 

(23) 

where {c q } k Zi are the positive coefficients associated to the 
Kq non-zero LOS atoms from the q- th source, and \c ql ^ 
are the positive coefficients associated to the I\ q i non-zero 
NLOS atoms between the source-sensor pair (q, l). Given R* 
is expressed as in (23) and given that (23) is an optimal atomic 
decomposition, i.e., its cost 

Q Kq Q Xj Kql 

c^EEtf’+EEEi'’ C4) 

q= 1 k—1 < 2=1 1=1 k =1 

is the smallest, then the set of locations for the g-th source 
associated with the optimal atomic decomposition is 

{p« for all k = 1,..., A',|, (25) 

the set of LOS propagation delays between source q and sensor 
l is 

{ T ; (p i k) ) : b ( q k \l) ± 0, for all fc = 1. K q } , (26) 

and the set of NLOS propagation delays between source q and 
sensor l is 

| T if } for all fc = 1,..., K q i |. (27) 

Next, a definition of correct recovery is provided. 

Definition 1. Given R* is expressed as in (23) and given 
that (23) is an optimal atomic decomposition, then the sources 
locations are correctly recovered if 

K q = 1 (28) 

= P q, (29) 

for q = 

Condition K q = 1 is required for all q because, obviously, 
it exists only one valid location for each source, and in such 
case p,' must match the true location of the <y-th source. 

In Table 1, the procedure for recovering the sources’ lo¬ 
cations from the received signals is summarized. In some 
specific cases, such as estimating frequencies from a mixture 
of complex sinusoids [26], some sophisticated techniques 
have been devised for minimizing the atomic norm, and then 
recover the optimal atomic decomposition of R*, thanks to the 
particular structure of the atomic set. However, in general, it is 
challenging to solve (22), because computing the atomic norm 
is not always straightforward. In Section V-B, an approximate 
method based on discretizing the atomic set is proposed for 


R 



Non-zero atoms < - > Proxy for locations 

Coefficients <- > Proxy for amplitudes 


Fig. 1. Flow diagram of the process for recovering the sources’ locations. 

simultaneously solving the atomic norm minimization problem 
(22) and recovering the optimal atomic decomposition (23). 
Before delving into the details on how to actually solve 
problem (22) and find the optimal atomic decomposition as 
expressed in (23), it is shown in the next section that tuning 
the parameters {u q }® =1 appropriately is critical to the correct 
recovery of the sources’ locations. 

A. Guarantee for Correct Recovery of the Sources ’ Locations 

In this part are developed guarantees for correct recovery of 
the sources’ location in the sense of Definition 1. To ensure an 
identifiable signal model, the following assumption is made: 

Assumption 1. For each sensor, signal model (5) is identi¬ 
fiable in the sense that the observed data is explained by a 
unique set of delays, T;(p) for q = 1,..., Q, and for 
q = 1,..., Q and m = 1,..., M q i. 

Further, to develop correct recovery guarantees, we assume 
noiseless observations, in which case the solution to (22) 
is trivially R* = R. However, as shown in the numerical 
section, the theoretical results obtained in this section are also 
meaningful in the presence of noise. The key properties that 
are exploited to obtain guarantees are: 

1) LOS signal paths associated with a source have a com¬ 
mon location (see (2)). 

2) NLOS signal paths are local to sensors (see (4)). 

To formalize the notion that LOS paths emitted by a source 
have a common location, we introduce the notion of location 
consistency: 

Definition 2. A location p is said to be consistent with X 
paths (LOS or NLOS), or vice-versa, if the propagation delays 
of such paths, say t\,..., tx, satisfy 

t x = T ix (p) for x = 1,..., X, (30) 

where {li,...,lx} Q {1, ...,L} are the indexes of the 
destination sensors of the X paths, and T; a . (p) is the delay 
of the direct path between location p and sensor l x . 

In order to find the sources’ locations exploiting the notion 
of consistency in Definition 2, the following assumptions are 
made. 
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Assumption 2. The number of LOS paths from source q, S q , 
is known. 

By its very nature, a source location cannot be consistent 
with any NLOS. Thus the location of the g-th source is 
consistent with exactly S q paths. 

Assumption 3. Only the true location of the q-th source is 
consistent with S q paths emitted by the q-th source. 


For the proof of Lemma 2, see Appendix B. The interpre¬ 
tation of this lemma is that given a solution that does not 
produce a location for the q-th source, and if condition (32) 
is met, there exists another lower cost solution that produces 
a location for the q-th source. 

The two lemmas lead directly to the following theorem 
establishing the guarantee for correct recovery of the sources’ 
locations. 


By Assumptions 2 and 3, given a source with a known 
emitted waveform and a known number S of LOS paths, its 
true location is the only one consistent with S paths. 

From (23) and Definition 1, the solution containing the 
true locations of the sources is associated with the optimal 
atomic decomposition. However, from (18) and the definition 
of atoms, namely, LOS atoms (15) and NLOS atoms (16), the 
optimal atomic decomposition is parameterized by the norm 
of the amplitudes in the LOS atoms u q (15). For given data 
R\ decreasing u q has to be balanced by an increase in the 
coefficients of the LOS atoms, thus raising their contribution 
to the cost C (24). Put another way, different values of u q lead 
to different explanations of the data R* manifested as different 
optimal atomic decompositions, and thus corresponding to 
different solutions of the source localization problem. We 
seek to determine which values of parameters u q ensure that 
the corresponding optimal atomic decomposition results in 
locations that are consistent with the number of paths indicated 
by Assumption 2. This in turn guarantees that these are the true 
sources’ locations. The next lemma establishes the condition 
on u q under which a location associated with the optimal 
atomic decomposition is also consistent with the number of 
LOS paths. 

Lemma 1. Given a known number of LOS paths S q of the 
q-th source, if parameter u q satisfies 

Uq < (31) 

then any location (for the q-th source) associated to the 
optimal atomic decomposition (25) is consistent, in the sense 
of Definition 2, with S q or more paths. 

For the proof of Lemma 1, see Appendix A. The interpre¬ 
tation of this lemma is that given a solution that produces 
a location with less than S q paths, and if condition (31) is 
met, there exists another lower cost solution, implying that a 
solution with fewer than S q paths cannot be optimal. 

The previous lemma guarantees that any location associated 
with the optimal atomic decomposition is consistent with S q 
paths. The next lemma establishes the condition on u q that 
ensures that at least one location is associated with the optimal 
atomic decomposition. 


Lemma 2. Given a known number of LOS paths S q of the 
q-th source, if parameter u q satisfies 


u q > 


1 

W*' 


(32) 


Theorem 1. A sufficient condition for the correct recovery in 
the sense of Definition 1 of the sources’ locations is that 


1 

7K 


< U q < 


1 

y/S^l 


(33) 


for all q. 


Proof: If u q > fiyGf for all q, by Lemma 2, at least one 
location is associated to the optimal atomic decomposition for 
each source. By Assumption 2, the number of LOS paths S q 
is known for each source q. Therefore, if u q is chosen such 
that u q < V sjs q -1 for all q, then by Lemma 1, the locations 
associated to the optimal atomic decomposition for the source 
q are consistent with S q or more paths. However, according 
to Assumption 3, only the location of the source is consistent 
with S q or more paths, thus completing the proof. ■ 

The interpretation of Theorem 1 is that the LOS atoms 
should be large enough to guarantee at least one LOS solution, 
but small enough to guarantee that the solution is the correct 
one. A numerical examples illustrates Theorem 1. Let the 
search area be of size 200 m x 200 m and centered around the 
origin of the coordinate system. A single source is positioned 
at (20 m, 30 m) and 5 sensors are positioned at coordinates 
(40 m, -40 m), (-40 m, -40 m), (-40 m, 40 m), (40 m, 40 m) 
and (0 m, 0 m). All sensors receive a LOS path except for 
the sensor located at (40 m, -40 m). Therefore, the number of 
LOS paths is Si = 4. In addition, the sensor at the origin 
receives a NLOS path whose path length is 91m. The goal is 
to compute the probability of correct recovery in the sense of 
Definition 1 as a function of u\ and under the conditions of 
Theorem 1, i.e., the noiseless case. The implementation of the 
procedure leading to Fig. 2 is discussed in Section V-B. To 
estimate the probability of correct recovery, the experiment 
is repeated 1000 times, and in each experiment the emitted 
waveform as well as the amplitudes of the LOS and NLOS 
paths are chosen randomly. The exact model for generating 
the waveforms, as well as other parameters is the same as 
the one detailed in Section VII. Figure 2 plots the probability 
of correct recovery versus parameter v which is defined as 
v = (V«i) 2 . Theorem 1 guarantees a correct solution if 


Si — 1 < v < S\. 


(34) 


As it can be seen in Fig. 2, for values of v within the interval 
]3,4[, the probability of correct recovery is close to one, 
whereas it is smaller for other values. 


B. Practical Implementation: Discretization of the Atomic Set 


then at least one location (for the q-th source) is associated Remind the reader that in general, when noise is present 
to the optimal atomic decomposition. the process for recovering the sources locations follows Fig. 1. 
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Fig. 2. Probability of correct recovery in the sense of Definition 1 in the 
absence of noise. 


The most straightforward method for solving problem (22) and 
obtain its optimal atomic decomposition (23) is to substitute 
the atomic norm in the objective function ( 22 a) by its defi¬ 
nition ( 21 ), and optimize over the set of positive coefficients 
{<;(*>}. However, such approach yields an infinite-dimensional 
problem because the number of atoms is infinite. Except for 
some particular cases, like recovering frequencies of mixtures 
of sinusoids [26, 27], it is in general very challenging to 
optimize infinite-dimensional convex problems. In [28], it is 
advocated that dictionaries whose atoms depend on continuous 
parameters are discretized. For instance, the NLOS atoms 
(16) depend on a delay, which is by definition within the 
interval [ 0 , r max ], and can be discretized into a grid of discrete 
delays such as (9). In [28], it is proven that the optimization 
problem based on the discretized atomic set converges to the 
original problem (22) as the grid finesse increases. Indeed, grid 
refinement approaches can be found in some signal processing 
applications such as delay estimation [23], direction-of-arrival 
estimation [29-31] or direct localization of sources [7], 

The atomic set is composed of LOS (15) and NLOS atoms 
(16). The LOS atoms are parametrized by the location of the 
source whereas the NLOS atoms are parametrized by their 
propagation delays. Therefore, two different types of grids 
need to be created: one grid of locations and one grid of delays. 
The propagation delays of the NLOS paths vary between 0 and 
77 nax- Upon discretizing the interval of delays with a resolution 
of T re s 

T> = { 0 , r res , ■ • • , 


(35) 


a new set of NLOS atoms is obtained 

Q L 


AnlOS = UU{ N q i (r) : 0 < (j) < 27r,r G pj. (36) 


q=l 1=1 


Similarly, upon discretizing the search area S into a uniform 
grid of squared cells whose center points are 


Q = {0!,- - ■, 0|e|}, (37) 


where the grid resolution is defined as d Tes = min ^ ||0j — 
0j || 2 , a new set of LOS atoms is obtained 


Aos= U {Mb,p) :beC L ,peScR 2 ,||b || 2 = u 9 }. 

9=1 

(38) 

The discrete atomic set including the LOS and NLOS atoms 
is 

-4 = Alos U -4nlos7 (39) 


and the atomic norm induced by A has the same expression 
than in ( 21 ) 


R 

_ = inf < 


A. c( fc )>0 j 


c (fc) : R = c W A (fc) , A (fe ) eA\, 


(40) 


except for the fact that A has been replaced by A. By 
expressing the generic atoms A 1 ' 1 '' 1 in (40) as LOS or NLOS 
atoms, the new atomic norm || • || j may be cast as 


R 


( Q |6| Q L \T>\ 

x= EE4 fl) +EEE4? (41a) 

c i >0 [ g=lfc=1 5=1 ( =1 9 =1 


with coefficients Cg 9) and such that 


+ 


Q 161 

R = EE4 5)l 9 (H 9) 7 0 fi 

9=1 9=1 

Q L \T>\ \ 

+EEE L (41b) 

q— 1 1=1 d=l ) 


and bq 9 -* may be any vector such that ||b ^ 9 ' ) ||2 = u q for all 
q and g. By replacing the LOS and NLOS atoms in (41b) by 
their definitions (15)-( 16), constraint (41b) becomes 


Q 161 

U = EE C 9 9)fc 9 9) ( Z ) S 9 ( T l (0 9 )) + 

9=1 ff =l 

Q L \D\ 

+ E E E cfe *' 1 S 4 (( d - 1 ) T ») 

q=l 1=1 d=l 

for / = 1,..., L (42) 

where R = [ri • • • r^] and b^ 9 ^ = [hq 9 ^(l) • • ■ b q 9 \l)] T . Next, 
substituting the atomic norm || ■ ||^ instead of || • ||_4 in problem 
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(22) with (41a) and (42) yields 

Q \S\ Q L \v\ 

EE4 9) +EEE 


Sd) 


r (g) -i , 

°q ’°ql q= 1 fc=l 

l|b<%=u„ 

0<^f<27T 


9=1 i=l g=l 


(43a) 


s.t. 


X] ll r « - f ill2 < £ 


(43b) 


1 = 1 

Q \G\ 

r i =EE4 9)6 « 9) (Os g Cn(0 g )) + 

9=1 5=1 
Q L |©| 

+E E E c ? e * 0< " s 9 - 1 ) Tres ) 

9=1 ;=1 d=l 

for l = 1,..., L. 


(43c) 


Problem (43) is not convex because of the bilinear forms, 
Cq 9 ^bq 9 \l) and c^e l ^n , appearing in constraint (43c). This 
can be easily remedied by the following variable changes 

ci 9) bi s) = y(s) 


q q J q 


from which it follows that 


.( 5 ) 13 ( 9 ) 




— C (S) y — 
L q a Q 


y(») 

J q 


_ „(<*) _ 
~ C ql - 


M)t 




(44a) 

(44b) 

(45a) 

(45b) 


Combining (44) and (45) with (43c) and (43a), respectively, 
results in the following optimization problem 


Q \S\ 

EE 


y[ 9) 


mm 

v (s) 

9 = 19=1 

4 Ho 


Q L \v\ 

EEE 

q= 1 1=1 d=l 


M)( 


Z"q'{l) 


(46a) 


s -t- Ell r '-^ll2 ^ e 


1=1 

Q \G\ 

*1 = EE4 9) (0 s 9 ( T l(Og)) + 

9=1 9=1 
Q \v\ 

+E E 4 d) ( Z K (( d - i)^) 

g—1 d= 1 

for l = 1,..., L, 


(46b) 


(46c) 


which is convex and finite-dimensional. Problem (46) is equiv¬ 
alent to the latent group Lasso problem [32], and specific 
algorithms for solving (46) exist in the literature [33]. More¬ 
over, the problem also falls into the class of second-order 
cone programs (SOCP), a subfamily of convex problems, for 
which efficient algorithms are available [34]. In our case, the 
SOCP type of algorithms resulted in the fastest computational 
times. The variable yq 9 \l) represents the amplitude of a LOS 
paths from source q to sensor l with delay 7 ~i(Q g ), whereas 
the variable Zq d \l) represents the amplitude of a NLOS path 
from source q to sensor l with delay (d— l)T res . Let {yq 9 '*} and 
{zq (l)} be the solutions to problem (46). Then, the location 


of the q-th source is the grid location Q g for which 11y^ 11 2 
is larger than zero. Intuitively speaking, minimizing the term 
Y^ ( q=i SE EE \ z q 1 ^0‘)\ i* 1 ^ objective function (46a) 
induces a sparse number of NLOS paths, whereas minimizing 
SgLi llyg 9 ^II 2 induces a sparse number of sources’ 
locations. 


C. Estimation of the Number of LOS Sensors 

According to Theorem 1, we must fix u q to a value that 
satisfies 1 /\fs q < u q < 1 /y / S q -i for each source q , where S q 
is the number of sensors receiving a LOS component from 
the q- th source. Hence, u q must be set to u q = 1 / y /S q -n for a 
parameter /i g]0,1[. For instance, it has been observed that a 
satisfactory choice is /z = 0.2 as it led to the best probability 
of correct recovery for all experiments in Section VII. In 
this section, we propose a method for estimating the sources 
locations that not only does not require a priori knowledge on 
the number of LOS sensors S q , but in fact estimates them. 
The method works as follows. We start by assuming that all 
sensors receive a LOS component from all sources, S q = L 
for all q, and set u q such that it satisfies (33). Then problem 
(46) is solved. According to Lemma 1, the sources’ locations 
associated to the optimal atomic decomposition for the q-th 
source are consistent with at least S q paths. However, by 
Assumption 3, no location is consistent with more than S q 
paths. Therefore, if the number of LOS sensors (S q > S q ) had 
been overestimated, no location would be obtained for source 
q. In the next step, S q is decreased by one for all those sources 
without a location estimate, and problem (46b) is solved again. 
These steps are repeated until a location is obtained for each 
source. The last value of S q is the estimated number LOS 
sensors for the q-th source. This method corresponds to steps 
10, 11, 20-28 of DLM’s algorithm described in Section VI. 


D. Spurious Locations 

It is observed in numerical simulations that when the 
sources are off-grid (p q (f Q for any q) and/or when the 
propagation delays of the paths are off-grid (r^ TO \ 77 ( p q ) f. V 
for any q, l, m ), then some spurious locations may be obtained 
from problem (46). This phenomenon is not new and it was 
studied in [35] in the case of delay estimation using the t\- 
norm (11). It was shown that if the propagation delay of a path 
is off-grid, a peak appears around such propagation delay but 
also secondary peaks of much weaker strength appear further 
apart. 

To eliminate spurious locations, it is set a simple threshold 
criterion. Let y q l! be, for all q and g, the solution to problem 
(46). Ideally, for each source q, ||yg 9 '*|| is zero for all q except 
if 0 g matches the location of the source. However, in practice, 
for a given source q, problem (46) may produce some spurious 
locations, in which case ||yg 9 '*j| may be different than zero 
for more than a single value of q. If 0, ; is the true location 
of the q-th source, then the S q largest components of y q 9 ^ 
are the amplitudes of the S q LOS paths. In contrast, if Q g 
is a spurious location, we have observed through numerical 
experimentation that some of the S q largest entries will be 
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approximately zero. Denote y g 9 ^ the vector with the same 
components than y q 9 \ but sorted in descending order, i.e., 
|^(s)4-(^)| > ... > We propose that, for a given 

source q, all locations whose S q strongest components do no 
satisfy 


(s)J- 


(s q ) 


> AT. 


(47) 


are dismissed. Here, S q is the number of LOS paths assumed 
for source q as explained in more detail in Section V-C, 
parameter A is the largest signal strength of a LOS or NLOS 
path 



and T is a value smaller than 1. For instance, in the simulations 
it was used T = 1 /.30, so that all locations whose signal 
strengths are 201og 10 (30) s=s 30dB weaker than the strongest 
path are discarded. If after the threshold criterion (47) one 
or more locations still remain for the q -th source, then the 
location with the largest strength is picked 


P q = ®9 


g = arg max 

9 




( 4 ) 


(49) 


It is important to not skip (47), and apply (49) directly. 
As explained in Section V-C, the proposed technique works 
by initially assuming that the number of LOS paths for 
the g-th source is S q = L and if no location is obtained, 
then successively decreasing S q until problem (46) outputs a 
location. However, if the threshold criterion (47) is skipped 
and S q > S q , a spurious location may be erroneously selected 
as the correct source location instead of concluding that there 
is no location and that S q needs to be decreased. 


E. Tuning Parameter e 

Parameter e in optimization problem (46) constraints the 
fitting error between the received signals and the estimated 
signals. Such a parameter is set so that the received signals 
without noise are a feasible solution. Let f; be the noiseless 
received signal at sensor l, then we require that 
L L 

ll r ' ~ Tilla = ^ e - (5 °) 

l=i i=i 

If e is chosen too small, then it can happen that J^i=i || w i ||2 ft 
e, thus excluding the noiseless signals from the set of possible 
solutions. Because the noise {w;}^ are random independent 
complex Gaussian vectors of length N, it follows that the 
error normalized by the noise variance 2 a~ 2 J2r= i Il w z|l 2 is 
a Chi-square random variable with 2 NL degrees of freedom. 
Thus, parameter e must be set to a large enough value so that 
EtilMI < e is satisfied with high probability, e.g., 

Pr (E llwilla < ^ = 7- (51) 


Let F(x, k) be the cumulative distribution function of the chi- 
squared distribution with k degrees of freedom evaluated at x 


and F its inverse function. Then 
_2 


^F" 1 (7,2JVL). 


At low signal-to-noise ratio (SNR), it is possible that the 
energy of the received signals is too low compared to the 
energy of the noise causing that IIdII! — e -1* 1 suc h case 

problem (46) has the trivial solution y q 9 \l) = z q 9 \l) = 0 
for all q, l, g and d, and it will not output any locations, 
if Ef=iM2 < e, we propose to estimate the locations by 
finding the LOS signals that have the highest correlation with 
the received signals: 


L 

p q = arg max 

1=1 


S q (n (p))l-;| 2 
\Sq{n ( p ))|| 2 


(53) 


F. Grid Refinement 

The computational complexity of minimizing the second- 
order cone problem (46) is 0((Q\Q\ + QL\U\) 3 - 5 ) [36]. To 
lower it we propose a recursive grid refinement procedure 
inspired by the ones in [29-31]. The optimization problem 
(46) employs a grid of delays in order to estimate the NLOS 
paths between every source-sensor pair, and a grid of locations 
in order to estimate the location of every source. In total, 
there are Q grids of location and QL grids of delays. In 
comparison to previous grid refinement approaches, ours is 
a more complex due to the two different types of grids used 
to explain the observed data. The idea behind a grid refinement 
procedure is to start with a coarse grid(s) and refine each 
grid only around the active points. Let T res and d res be the 
grid resolutions we wish to achieve in the grids of delays and 
locations, respectively, and suppose that in order to lower the 
computational complexity, the grids are refined R times. If the 
resolution of the grids is increased by a factor of two at every 
step, then the grids resolutions at each step are 

7-res,r = 2 fl-r T res for r = 1,..., R (54) 

fires,r = 2 R ~ r d les for v = 1,..., R. (55) 

Let V q i tr be the grid of delays for the source-sensor pair (q, l) 
at step r, and Q qr the grid of locations for source q. At the 
first step (r = 1), the continuous set of delays [0, r max ] is 
discretized with resolution T res ,i 

77g/,i = {fT re s,i € [0, T max ] : i € , (56a) 

and the search area S is discretized uniformly with resolution 

fires, 1 

Qq,i = (jj e S : i, j e zj , (56b) 

where Z is the set of integers. Consider step r, and let the 
active propagation delays between the source-sensor pair (q, l) 
be : m = 1,..., M q ^ r }, and the active locations for 

source q be {pq'r : m = 1,..., K q r }. Then, the grids at step 
r + 1 include the previous active delays and locations plus 
some neighbor points. For instance, in addition to the active 
delays and locations, we include two points at the left and 
right of the active delays 

Mql,r 

v qhr+1 = [J {t^+ iT res ,r+1 : * =-2,-1,0,1,2|, 

m—1 

(57a) 


e = 


(52) 
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Step 1 


Step 2 


Step 3 


I . .. 


N 


Estimate: r* 


r* 4- t '' 

T f T r( 


Estimate: 0* 


T 

I 

r~ 

. 1 .... 

i 

I--T 

I 

-*l _J_ _ 

e* + 


i 

i 


d (2) 

“res 

-d (2) 

“res. 


-*r _ 

i 






0 * + 


d (2) ‘ 
“res 

-d (2) 

“res. 


d (3) l 

“res 


Fig. 3. Illustration of three steps of a grid refinement procedure. The top 
image shows the grid refinement for the delays between a hypothetical source 
and sensor, and the bottom image shows the grid refinement for the locations 
of some hypothetical source, for r = 1, 2, 3. The dots point out the position 
of a non-zero delay and location as a result of optimizing problem (58). The 
positions of such non-zeros are progressively refined at each step. 


and all points within distance 2d res ,r+i in the x- or y-axis of 
the active locations 


Q q ,r+ 1= U (p^+tW+i (•) :*,j = -2,-l ) 0,l,2. 

(57b) 

For a more intuitive picture on the grid refinement procedure 
see the examples in Figure 3 with three steps. Next, problem 
(46) is solved again but only for the new grid points: 


mm 

WT 

W d) } 


Q 

E E 

0—1 9- 

9g£Gq,r +1 


Yq 9) 


Q L 

-+EE E 

q= 1 1=1 d: 

(d— l)r„ s 


4 d) (0 


s.t. 


£ 


1=1 


(58a) 

(58b) 


r*=E E y q 9) ( l ) s q( r i( Q g)) 


q= 1 a- 

0 s e0g,r+1 

Q 


E E 4 d) ( l ) s q i( d ^ ^ 


9=1 d: 

(d- 1)t« s 

£'Dql,r+l 


for l = 1 ,..., L. 


(58c) 


The process of refining the grids and solving problem (58) 
is repeated for the R steps. The proposed grid refinement 
procedure corresponds to steps 12-18 in DLM’s algorithm 
described in Section VI. 

In regards to the resolutions of the grids, instead of choosing 
the resolution of both types of grids completely independently. 


they are set according to 


CTre 


= d n 


for any r, 


(59) 


where c is the speed of light. 


VI. Algorithm 

In this section, it is presented the proposed DLM algorithm 
for source localization in multipath. The inputs to the DLM 
algorithm are the received signals {r /}^ =1 and the noise 
variance o' A w . The number of sensors L, sources Q and samples 
per sensor N are assumed known. The outputs of the algorithm 
are the source locations estimates {p, ; \, l= \ • The summary of 
the proposed algorithm for direct localization of RF sources 
in the presence of multipath is as follows: 

Input: L, Q, N, and <r 2 W . 

Parameters that need to be selected: S, r max , d res , T, 7 . 

Output: The source locations estimates {p q\® = i 

Procedure: 

1 : for sensor l where l = 1,..., L do 
2 : Estimate multipath TOA’s {f^} using [23] or any 

other delay estimation technique of choice. 

Estimate multipath amplitudes {5^} through (13). 
Reduce NLOS interference on the received signal r; 
through (14). 
end for 

Compute parameter e through (52). 

if Ei-i INI2 > e then 

Compute the initial coarse grids with (56) and (59). 
Initialize S q = L for q = 1,..., Q 
while p q = 0 for any q £ { 1 ,..., Q} do 


3: 

4: 

5: 

6 

7 

8 
9 

10 

11 

12 

13: 

14: 

15: 

16 

17 

18 

19 

20 
21 

22 : 

23: 

24 

25 

26 

27 

28 

29 

30 

31 

32 


for q = 1,... ,Q 


for r = 1,..., R do 

Optimize problem (58). Output: {y^, 1 } and 

{^( 01 - 

Find the active locations {j>q'r } and the active 
delays {t^}. 

if r / 1 then 

Refine the grid with (57) and (59). 

end if 
end for 

Compute A through (48). 
for q = 1,..., Q do 

if any locations are active for the q-th source 
and such locations satisfy (47) then 

Estimate the location of the ( 7 -th source 

through (49). 

else if S q > 1 then 

S q <- S q - 1 

else 

Estimate the location of the ( 7 -th source 

through (53). 

end if 
end for 
end while 
else 

Recover sources’ locations through (53). 

end if 
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Fig. 4. Map with the locations of the sensors and source used in many of 
the experiments in Section VII. 


VII. Numerical Results 

In this section, we illustrate the performance of the lo¬ 
calization method by numerical examples, and compare it to 
other existing techniques via Monte Carlo simulations. In all 
examples, the sources and sensors are positioned within a 
square area of 200 m x 200 m, which is divided into a grid 
of 1 m x 1 m cells, thus resulting in 40,000 cells. Unless 
stated otherwise, we simulate a scenario containing one source 
positioned at coordinates (20 m,30 m) and 5 sensors positioned 
at coordinates (40 m, -55 m), (-45 m, -40 m), (-50 m, 55 m), 
(60 m, 60 m) and (5 m, 0 m) as pictured in Fig. 4. The signals 
emitted by the sources are drawn from a white Gaussian pro¬ 
cess and filtered so that their passband bandwidth is 10 MHz. If 
multiple sources, such as in the experiment of Section VII-F, 
the waveforms are generated independently, thus the cross¬ 
correlation between signals from different sources is low but 
not necessarily zero. All sensors are time-synchronized and 
sample the received signals at a 20 MHz frequency for a total 
time of 5 ps, thus each sensor observes 100 samples. For each 
source, we define the SNR per observation time as 


components at all sensors are modelled by a Poisson process. 
The mean inter-arrival time is set to 0.2 ps, and the average 
power P of a NLOS arrival at sensor l is governed by the 
power delay profile (PDP) 

Pi{t) = exp I —-—-— I (61) 

\ £rms J 

where t is the arrival time of the NLOS component, is 
the arrival time of the LOS path and t rms is the root mean 
square (rms) delay spread. An exponential PDP assigns smaller 
power to later arrivals. Unless otherwise stated, all LOS paths 
have normalized unit power. In multipath environments, it is 
possible that some sensors have their LOS blocked, thus at 
each Monte Carlo repetition one randomly selected sensor 
among the five receives no LOS component. 

The figures compare the performance of the following two 
direct localization techniques: 

1) DLM — The proposed technique. 

2) DPD — Direct Position Determination as originally 
propose in [6] for AWGN channels. 

3) DPD with NLOS mitigation — In this variation, DPD 
is preceded by the NLOS mitigation method introduced 
in Section IV. The goal is to show that DLM outper¬ 
forms this variation of DPD, to demonstrate that DLM’s 
high accuracy is not due only to such NLOS interference 
mitigation method. 

4) Indirect, CS TOA — Indirect localization comprises a 
two-step process. In a first step, the TOA of the first path 
at each sensor is estimated by a delay estimation method 
based on compressive sensing (CS) [23]; in a second 
step, multilateration is performed using the well-known 
method developed by Chen [3] to mitigate the problem 
of potential LOS blockage on sensors. 

5) Indirect, matched filter TOA — Same as previous 
indirect technique, except that TOA’s are estimated by a 
threshold-based matched filter. 

To solve the conic problem in DLM (step 13 of DLM’s 
algorithm described in Section 13) and in CS TOA, we utilize 
the Mosek solver [38]. The bandwidth of the emitted signals 
limits the localization accuracy, and it is known that the 
ranging resolution is approximately 


/ NLP los \ 

SNR = 10 log 10 (-^ , (60) 

V / 

where N is the number of observations per sensor, L is the 
number of sensors, Plos is the power of a LOS component, 
and is the variance of the sampled noise. According to 
[37], in urban and suburban areas, the signal strengths of 
LOS and NLOS paths may be modeled as random variables 
with log-normal distribution. It follows that the channel tap 
powers expressed in dB are random variables with normal 
distribution. For our simulations, we set the standard deviation 
of the tap powers to lOdB. All multipath experiments simulate 
Turin’s urban channel model [37]. The arrival times of NLOS 


where c is the speed of light and B is the signal bandwidth. 
For the particular case of a 10 MHz bandwidth, the waveform 
ranging resolution is then 30 m. Also, we define the probability 
of correct recovery for the case of a single source as 

1 Z 

p c = ^E 1 ( |p_ p (z) i <c )’ (63) 

where p is the true source’s location, Z is the number of times 
that the experiment is repeated, p (z ' ) is the source’s location 
estimate for the z-th repetition, and I (■) is the indicator 
function. Unless otherwise stated, the error is set to ( = r /z, 
which is a value smaller than the ranging resolution r. In some 
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Fig. 5. Root mean square error vs. SNR for the scenario in Fig. 4 when no 
multipath is present. 



Fig. 6. Probability of correct recovery vs. SNR for the scenario in Fig. 4 
when no multipath is present. 



Fig. 7. Root mean square error vs. SNR for the scenario in Fig. 4 in a 
multipath environment. 



Fig. 8. Probability of correct recovery vs. SNR for the scenario in Fig. 4 in 
a multipath environment. 


of the tests, it is plotted the normalized root mean square error 
1 


rMSE = 


A 


1 


J2(p-p (z) ) ■ 


(64) 


Z=1 


All experiments are repeated 1000 times, i.e., Z = 1000. 


probability of recovery because essentially both techniques, in 
the absence of multipath, look up for the location whose LOS 
signals correlate the most with the received signals. DPD and 
DLM perform substantially better in comparison to indirect 
techniques as it is expected from the theory. 


A. Performance in the Absence of Multipath 

This experiment’s purpose is to validate that DLM per¬ 
forms optimally in the absence of multipath, i.e., its accuracy 
matches that of the DPD, which was shown to be optimal (see 
[6]). All five sensors receive LOS components, and Turin’s 
channel model does not apply here, since there are no NLOS 
paths. Ligures 5 and 6 plot the rMSE and the probability 
of correct recovery, respectively. DPD and DPD with NLOS 
mitigation are plotted together because their performance is 
exactly the same in the absence of multipath. As it can be ob¬ 
served, DPD and DLM perfom equally in terms of rMSE and 


B. Performance in Multipath 

In this example is simulated the multipath channel model 
described at the top of this section. The rMSE and the 
probability of correct recovery vs. SNR are plotted in Ligs. 7 
and 8, respectively. Observe in Ligs. 7 and 8 that DPD fails to 
localize the sources irrespective of the SNR due to the fact that 
it is not designed for multipath. Also, the indirect technique 
relying on estimating by matched filter the TOA of the first 
arrival, does not perform much better than DPD because 
matched filter suffers from severe bias when multiple arrivals 
overlap in time. Interestingly, it seems as if DLM does not 
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Fig. 9. Probability of correct recovery vs. error for the scenario in Fig. 4 Fig. 10. Probability of correct recovery vs. rms delay spread for the scenario 
for a 30 dB SNR. in Fig. 4 in a multipath environment for a 30 dB SNR. 


perform better, in terms of rMSE, than the indirect technique 
employing CS TOA estimates. In Fig. 9, the probability of 
correct recovery (63) is plotted for different errors ranging 
from 0 to 2 r for an SNR value of 30dB. DLM achieves a 
high probability of correct recovery for much smaller errors 
than the other methods. For instance, DFM’s probability of 
correct recovery is 0.9 for an error smaller than 0.4r, whereas 
for the indirect technique with CS TOA, such probability is 
only achieved when the error is 0.9r. The other techniques 
perform substantially worse than DFM, and in fact, they never 
achieve a probability of recovery close to one even when very 
large errors are allowed. In summary, DFM can achieve a 
high probability of recovery for very small errors. In terms of 
rMSE, DFM and the indirect technique employing CS TOA 
estimates perform similarly, because in the rMSE metric small 
errors have a much smaller impact compared to the large 
errors. Hence, in the next experiments, we focus only on the 
probability of correct recovery. 

C. Probability of Correct Recovery vs. Delay Spread 

The considered channel model depends on the rms delay 
spread, which determines the interval between the EOS com¬ 
ponent and the last arriving NFOS component. In general, 
larger delay spreads imply more multipath that make the 
localization more challenging. In Fig. 10, the probability of 
correct recovery is plotted for an rms delay spread ranging 
from 0 to 0.6 ps at 30 dB SNR. At high-SNR and at a zero 
delay spread all localization techniques perform similarly. 
However, as soon as the rms delay spread increases by 
a little as 0.2 ps, DPD’s performance drops markedly. The 
techniques specifically designed for multipath channels, such 
as the indirect technique based on CS TOA estimates and 
DFM, degrade very slightly as the rms delay spread increases. 
DFM outperforms all other techniques and is capable of 
recovering the sources locations with a high probability of 
correct recovery irrespective of the delay spread. 



Fig. 11. The left axis plots the probability of correct recovery and the right 
axis the mean elapsed time for running DLM’s Stage 2, vs. the number of 
grid refinement steps. The SNR is fixed at 30 dB. 

D. Probability of Correct Recovery vs. Number of Grid Re¬ 
finement Steps 

The purpose of the grid refinement procedure introduced 
in Section V-F is to reduce the computational complexity of 
DFM, while maintaining the localization accuracy. Figure 11 
plots the probability of correct recovery (square marker) and 
the DFM’s mean elapsed time at Stage 2 (circle marker), 
versus the number of grid refinement steps. The SNR is fixed at 
30 dB. DFM is run on a computer with an Intel Xeon processor 
at 2.8 GHz with 4 GB of RAM memory. Perhaps surprisingly, 
the probability of correct recovery remains almost constant 
irrespective of the number of steps. The lowest computational 
time is 5 s and is obtained for five grid refinement steps. The 
number of grid steps that results in the lowest computational 
time depends on many factors such as number of grid points, 
efficiency of the conic solver, particular scenario and so forth. 
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Fig. 12. Probability of correct recovery vs. the number of LOS sensors for Fig. 13. Average probability of correct recovery vs. the number of sources 
a 30 dB SNR. for a 30 dB SNR. 


Thus, in general, the optimum number of steps must be found 
by in situ testing. 


E. Number of LOS sensors 

The information about the sources’ locations is carried on 
the LOS components (see signal model (5)). This experiment 
evaluates DLM’s probability of correct recovery versus the 
number of sensors receiving a LOS path. Since the setup of 
Fig. 4 includes five sensors, the number of LOS sensors is 
varied between one and five. As expected. Fig. 12 shows that 
a larger number of LOS sensors results in better localization 
accuracy. For the cases where there is only one or two LOS 
sensor, the probability of correct recovery drops drastically 
because, in general, the minimum number of LOS sensors 
required for unambiguous TOA-based localization is three. 


F. Multiple Sources 

In this example is evaluated the probability of correct recov¬ 
ery of multiple sources emitting different signals overlapping 
in the time and frequency domain. The SNR is fixed at 30 dB. 
The definition of the probability of correct recovery defined in 
(63) was for a single source. In the case of multiple sources, 
we define the average probability of correct recovery 


Pav 


1 

ZQ 


Z Q 



(65) 


~ (z) 

where p q is the true location of the ry-th source, pf is its 
estimate, and £ is the error set to £ = r /3 where r is the 
waveform’s ranging resolution as defined in (62). In Fig. 13, 
it is shown how the average probability of correct recovery 
degrades as the number of sources increases. This is expected 
because the signals from different sources interfere with each 
other. Nonetheless, we can observe that DLM outperforms all 
other localization techniques when localizing multiple sources. 


VIII. Conclusions 

By combining concepts from compressive sensing and direct 
localization, we have developed a novel direct localization 
technique for mutliple sources in the presence of multipath 
(DLM). This technique assumes the emitted waveforms are 
known but requires no prior information on the channel. In 
fact, its localization accuracy is almost constant irrespective of 
the delay spread of the channel. At the core of our technique 
lies an optimization problem that recovers the locations of the 
sources with high accuracy by exploiting properties that are 
different for LOS and NLOS paths. It is shown theoretically 
how to set the algorithm’s parameters to guarantee successful 
recovery including a parameter that determines the relative 
contributions of the LOS and NLOS components to the cost 
function. Contrary to indirect techniques, the proposed tech¬ 
nique is capable of localizing sources with an accuracy beyond 
that of the signal resolution, with high probability. In absence 
of multipath, DLM’s accuracy matches that of the maximum 
likelihood estimator of the sources’ locations. In the presence 
of multipath, DLM’s accuracy is better than the accuracy of 
all compared methods and can find the sources’ location even 
when some sensors suffer from LOS blockage. The gain in 
localization accuracy does not come for free, as DLM requires 
larger computational resources than previous techniques. To 
this end, we propose a grid refinement procedure which sub¬ 
stantially reduces the computational complexity. Nonetheless, 
this should be less of a burden as computational power keeps 
increasing and second-order cone program solvers become 
more efficient. DLM’s high accuracy is validated by extensive 
numerical simulations. 


Appendix A 
Proof of Lemma 1 

Let an atomic decomposition of R be (23). The goal of the 
proof is to show that all locations, p!/' ; for q = 1 ,.... Q and 
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k = 1,..., K g are consistent with S q or more paths if 



From (23), the signal at the Z-th sensor is 

r< = E E c^6W(Z) s,(r,(pf)) + 

<7=1 /c—1 

Q ^qZ 

=,(’■?)■ (67) 

(?=1 fc=l 

By Assumption 3, T;(p|/' : ' ) ) is a true propagation if b q k \l) ^ 0. 
Therefore, if bj/^ has S q or more non-zero entries, according 
to Definition 2, p q is consistent with S q or more paths. It is 
left to prove that ||bq fe '*||o > S q . The proof is by contradiction. 
For instance, 

b^ 1} < Si. (68) 

o 

and let the atomic decomposition (23) in which the atom 
Li is replaced by ||bq^||o NLOS atoms as fol¬ 

lows 

L 1 (b( 1) ,p^)= £ 16^(01 N 1X (ti (pi 1} )) ■ (69) 

1 = 1 

Consider now the two decompositions (23) and the one 
obtained with (69). The costs of the two decompositions differ 
only in the coefficients of the atoms shown in (69). Ignoring 
the common atoms, the cost of decomposition (23) is c^\ 
whereas the cost of decomposition obtained from combining 
(69) with (23) is 

<4 1} E K } (D|. (70) 

1 = 1 

0 

Normalizing the two costs by cj 1 , and if (23), which by ( 68 ) 
has a location pj 1 ' 1 with less than S q paths, is optimal, then 

L 

1< \b ( i\l)\. (71) 

1=1 

* 4 1 } ( 0 # 0 

We show next that inequality (71) cannot be satisfied if 
||b ^ 1 ^||2 satisfies ( 66 ). Define the vector function l(b^) 
whose Z-th entry is one if 6^(0 7 ^ 0- an d 0 otherwise, and 
denote | • | the element-wise absolute value. Then the right 
hand side of (71) is 

£ | 6 A<)| = [l(b‘ ,, )] T |b< 1 >|. (72) 

1=1 

0 

and by the Cauchy-Schwarz inequality 



However, | b i 1} ||2 = u i, and by equation ( 66 ), ||bi 1) ||2 < 
•/v'S; i. Moreover, by assumption ( 68 ), ||b^||o < S\ — 1. 
Therefore, it follows 



which combined with (72) and (73) results in 

L 

E < !’ (75) 

1 = 1 

b?\l)* 0 

which contradicts (71). 

Appendix B 
Proof of Lemma 2 

Let (23) be an atomic decomposition of R. Recall that 
parameter K q is the number of locations associated to the 
optimal atomic decomposition for the q-th source. We aim to 
prove that if parameter u q 

K ) IL =u « > ^’ (76) 

then the optimal decomposition has K q > 1. The proof is by 
contradiction. Let K\ = 0, then a presumed optimal atomic 
decomposition (23) simplifies to 

»*EE 4‘ A (H‘>. pf>) +EEE (k?) ■ 

q= 2 k= 1 q= 1 1=1 k= 1 

(77) 

From (77), the signal at the Z-th sensor is 

d = E E 4 fc) ^ fc) (ZK^(p^)) + 

9=2 fc=1 

+££ 4 ! 7 ^ s \(}«')■ <™ 

q—1 k= 1 

Notice that the first summation begins with q = 2, because 
K\ = 0. By Assumption 3, 



are the true propagation delays of the paths between source 1 
and sensor Z. By Assumption 2, there are S\ LOS paths from 
source 1. Let 

{Zi,-.-,Zs 1 }C{l,...,L} (80) 

be the indexes of the destination sensors of such LOS paths, 
and let r ^ 1 in (79) be the propagation delay corresponding to 
the LOS path between source 1 and sensor Z, i.e., 

T ii )=T ‘( Pi) for Z € {Zi,...,Z 5l }. (81) 

We show next that there exists a decomposition different 
than (77) for which K\ > 1 and whose cost is lower, thus 
contradicting the assumption that (77) is optimal. According 
to (15) and (16), the sum of NLOS atoms with delays 
t -3 '' for l £ {li,... ,lg 1 j in the presumed optimal atomic 
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decomposition (77), i.e., i s } c u -^ii( T iV J )> can t> e 

expressed for any parameter c as 




E 

ZG{/i } 




KO 


ysic 

Ui 


Li (b,pi) 


+ £ ( c if — c ) Nil ( r if) 

iefii.-.-jisj} 


where b is 

m = 

Let C = C m i n 


A 1 ) 


_ j for l G {/i, 

0 otherwise. 


>Yi} 


defined by 

Cmin = mm 


„( 1 ) 


(82) 


(83) 


(84) 


ie{L,---7si} 

Next it is shown that the cost of the decomposition obtained 
by combining (82)-(84) with (77) is lower than the cost of the 
decomposition (77), contradicting the assumption that (77) is 
optimal. Notice the former decomposition includes the LOS 
atom Li(b,pi). The costs of the two decompositions differ 
only in the coefficients of the atoms shown in (82). Ignoring 
the common atoms, the cost of decomposition (77) is 


E 


„(i) 


(85) 

ze{L,-..4si} 

whereas the cost of the decomposition obtained from (82)-(84) 
is 

Vs~i Cmin 


Ml 


£ (■ 


J 1 ) - r 


( 86 ) 


Since (77) is presumed optimal, it means it must satisfy 


E 


(1) . yj ^l^min 


c \i — 


Ml 


£ (< 


.(i) 


(87) 

which after simplification leads to ui < 1 /y r Si, contradicting 
(76). 
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